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HIGHER-ORDER COMPACT SCHEMES FOR NUMERICAL 
SIMULATION OF INCOMPRESSIBLE FLOWS* 

ROBERT V. WILSON*, AYODEJI O. DEMl'REN* and MARK CARPENTER* 

Abstract. A higher order accurate numerical procedure has l>een developed for solving incompressible 
Navier-Stokes equations for 2D or 3D fluid flow problems. It is based on low-storage Runge-Kutta schemes 
for temporal discretization and fourth and sixth order compact finite-difference schemes for spatial discret- 
ization. The particular difficulty of satisfying the divergence-free velocity field required in incompressible 
fluid flow is resolved by solving a Poisson equation for pressure. It is demonstrated that for consistent glo- 
bal accuracy, it is necessary to employ the same order of accuracy in the discretization of the Poisson equa.- 
tion. Special care is also required to achieve the formal temporal accuracy of the Runge-Kutta schemes. 
The accuracy of the present procedure is demonstrated by application to several pertinent benchmark 
problems. 

Key words, compact schemes, incompressible flow simulation 

Subject classification. Fluid Mechanics 

1. Introduction. For direct numerical simulation (DNS) of fluid flow problems, it is generally 
accepted that higher-order accurate methods must be used ter minimize dissipation and dispersion errors. 
As the flow Reynolds number increases so do the ranges of temporal and spatial scales which must be 
resolved. Thus, the number of grid points- per- wavelength (PPWj recpiired by the numerical scheme for 
approximation of the flow equations to within acceptable tolerances of dissipation and dispersion errors 
effectively limits the smallest scales that can be computed accurately, and thereby also the maximum Rey- 
nolds number. Spectral methods require the fewest PPW, namely two, and are therefore ideal for computa- 
tions of flows with periodic boundary conditions. For more general flow problems finite-difference methods 
are desirable. Lele [1] has analyzed the resolution qualities of several finite-difference schemes. In general, 
resolution increased, i.e., fewer PPW, the larger the computational grid stencil, and implicit compact 
schemes had better resolution than regular explicit schemes of the same order of accuracy and computa- 
tional stencil. Further, high resolution properties could be improved by optimization of coefficients, but at 
the expense of the formal order of accuracy. A fourth-order compact scheme was devised with nearly the 
resolution quality of spectral methods, but retaining the flexibility of finite-difference methods. Haras and 
Ta’asan [2] analyzed the accuracy of various compact schemes, and demonstrated that resolution efficiency 
was not synonymous with minimum truncation error over the whole range of wavelengths present in a 
problem. The suggested that a scheme should be optimized for global accuracy rather than resolution effi- 
ciency. Hu et al. [3] have also shown that, in applications of interest to computational acoustics, temporal 
resolution could be improved by optimization of the coefficients of any multi-stage, Runge-Kutta, time- 
advancing approximation scheme. 
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NASI- 19480 while the first two authors were in residence at the Institute for Computer Applications in Science and Engineer- 
ing (1CASE), NASA Langley Research Center, Hampton, VA 23681-2199. Additional support was provided by the NASA 
Graduate St udent Research Program. 

. ^Department of Mechanical Engineering, Old Dominion University, Norfolk, Virginia 23529. 
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Iii large eddy simulation (LES) of turbulent flows, the goal is not to resolve all the scales in the flow, 
but only the larger scales. Effects of the unresolved smaller scales are approximated with sub-grid scale 
(SGS) models. Therefore, second-order, central-difference schemes are often used [4,5], A further justifica- 
tion is that effects of truncation errors may be comparable to uncertainties inherent in the SGS models. 
However, questions remain as to how large the "large” scales are, and how many PPW are required to 
resolve the smallest scale in the range. Further, what are the consequences of inadequate resolution of such 
second-order schemes. It is obvious that if LES is to be used in computational acoustics, dispersion errors 
are unacceptable [6]. Higher-order-accurate (greater than second-order) methods guarantee much better 
convergence towards grid independence, along with better wave-number resolution. In addition, implicit 
(compact) finite-difference schemes require narrower computational grid stencils, have better fine-scale res- 
olution and yield better global accuracy than explicit finite-difference schemes with the same formal order 
of accuracy. Therefore, the present study focuses on the use of liigher-order (fourth and sixth) compact 
schemes for the simulation of incompressible fluid flow problems. 

The lack of an evolution equation for the pressure presents particular difficulty in the computation of 
incompressible flows, which is absent in compressible flows. In the latter, the Navier-Stokes equations, 
along with the continuity equation and the energy equation present evolution equations for the five vari- 
ables, namely three momentum components, density and enthalpy, which can be advanced numerically 
with the same temporal and spatial discretization schemes. The pressure can then be obtained from an 
equation of state [7]. But in the former, an auxiliary equation has to be derived for the pressure which is 
then solved to satisfy the divergence-free velocity- field condition required for incompressibility [8,9,10]. Bell 
et al. [8] proposed a second-order projection method. Henshaw [9] derived a Poisson equation for pressure 
which was solved by a fourth-order explicit finite-difference method. It was found necessary to introduce a 
damping term to reduce divergence errors. In their finite-difference formulation for 2D problems, Joslin et 
al. [10] user! an influence matrix method to solve the Poisson equation for pressure. But this method is of 
order N for 2D problems and order \' :i for 3D problems, so that memory requirements and computation 
work quickly become prohibitive for 3D computations. If any of the directions has periodic boundaries, 
then mixed fmite-difference/spectral methods could be used [10,11], resulting in an Helmholtz equation for 
pressure which could be solved more cheaply. In the present study, interest is mainly in a flexible numerical 
scheme which could be used for non-regular computational domains and non-uniform computational grids, 
and with consistent treatment of all equations. Therefore, the Poisson equation for the pressure will be 
approximated with the same compact finite-difference scheme as used in the Navier-Stokes equations. 
Extension to irregular grids in physical space simply requires transformation of the equations onto a regu- 
lar grid in computational space. The metrics of the transformation must be computed with the same com- 
pact finite-difference scheme to guarantee a consistent level of accuracy. 


2. Analysis. 

2.1. Governing Differential Equations. The partial differential equations governing the incom- 
pressible fluid flow are the Navier-Stokes equations which can be written in Cartesian tensor form, for 
dimensionless variables as: 
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( 1 ) 


a«, 


du • 


-=: — h a — 

5/ J d.r 


dP t 7 d 2 », 

3r / Rc L dr ., dx .i 


where, u t are the Cartesian velocity components in the Cartesian coordinate directions x n Pis the pressure 
and Rr L is the Reynolds number based on the characteristic length, L. These equations must be solved in 
conjunction with the continuity equation: 


( 2 ) 


du 

i 

J7~. 

i 


= 0 


which expresses the divergence-free velocity condition. In 2D, i or j ~ 1,2 and in 3D, i or j = 1,2,3. Ein- 
stein’s summation rule for repeated indices is presumed. 

2.2. Temporal Discretization. The time advancement scheme for the momentum equations should 
possess several qualities such as low dispersion and dissipation errors over a wide range of step sizes, low 
memory storage requirements, and a relatively large stability envelope. A family of low-storage, RK 
schemes proposed by Williamson [12] possesses these desirable qualities. The schemes are low-storage in the 
sense that only two storage locations (one for the time derivative' and one for the variable' itself) are 
required for the time advancement. In comparison, a third-order fully implicit scheme would require four 
storage locations. 

The Navier-Stokes equations (1) are discretized temporally with explicit Runge-Kutta (RK) schemes, 
and spatially with implicit compact finite difference schemes. The momentum equation is advanced from 
time level, n, to n+1, in Q substages using either a third- or fourth-order explicit RK scheme. The 
advancement from substage, M , to M + 1 , is defined by: 


(3) 


M+l M 

y = U + 


b M+ 1 At^Hf 1 -8 X d 


M 


where At is the time step, £r 1+A is a coefficient of the RK scheme, and u- represents the Xj velocity com- 
ponent at the A/ t i substage. The substage, M = 0, is equivalent to the n t ^ x time level and the substage, M 

M M l c M M „Al- 


Q- 1, is equivalent to the nTl^ time level. The term, H . 


x 1 X 

= - M O U ■ + 7 ; O 

J X- l Re 


u . + a 

xx ■ i 


H. 


is the sum of the convection and diffusion terms, plus accumulation frail the previous sAib-stage. 5^ and 

8 , are compact first and second derivative operators, respectively, to be addressed in the next suf)-sec- 

tion.^The terms on the right hand side (RHS) of Eq. (3) are assumed known from the previous sub-stage or 

from the initial conditions at t = 0. The calculation of the pressure is accomplished by solving a Poisson 

equation at each sub-stage such that the continuity is enforced. Since the pressure, is calculated before 

A/ + I 

the advancement of Eq. (3), u i , can be calculated explicitly. 

The low-storage requirement is accomplished by continuously overwriting the storage location for the 
time derivatives and unknown variables at each sub-stage: 


(4) 


J.\/ M u U-t 
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( 5 ) 


AI+ 1 M , A/ + / A tt M 
u- <r- u- + b AtH- 


where fl l = fif 1 - [dP ^ 1 /dx-] and the notation <— is used to indicate that the storage locations, 
ft t , uf* are overwritten by, , uf** 1 , respectively. Tables 1 and 2 give values of the coefficients, 
a M and b M for some low-storage, three-stage- third-order and five-stage- fourth-order schemes, respectively. 


Table 1. Coefficients of a three-stage-third-order Runge-Kutta scheme, from Lowery and Reynolds [13] 


M 

a M 

b M 

1 

0 

0.500 

2 

-0.68301270 

0.91068360 

3 

-1.33333333 

0.36602540 


Table 2. Coefficients of a five-stage- fourth-order Runge-Kutta scheme, from Carpenter and Kennedy [14] 


M 

a M 

b M 

1 

0 

0.14965902 

2 

-0.41789047 

0.37921031 

3 

-1.19215169 

0.82295502 

4 

-1.69778469 

0.69945045 

5 

-1.51418344 

0.15305724 


The stability characteristics of RK schemes can be analyzed by considering the model equation: 

( 6 ) ^ t ) 

where $ is the generic unknown to be advanced in time and fl is the time derivative which contains the 
spatial terms of the governing equation. Equation (6) is transformed from physical space to wavenumber 
space by decomposing <|> into Fourier modes: 

(7) <|> = me ikx 

where <j>C) is the Fourier coefficient of <h , i = , and k is the wavenumber. Substituting Eq. (7) into 

Eq. (6) yields: 

(8) ^j<t> = A.<j> 
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where X (a complex number) is the Fourier symbol of the spatial operator fl . The RK scheme is used to 
expand the term on the LHS of Eq. (8). This gives the amplification factor, G = <(> / for the third- 

order scheme as: 


( 9 ) 


G = 1 + (XAt) + ^(XAt) 2 + t(KAt) 3 


It can he shown that all three-stage, third-order RK schemes have the same amplification factor given 
in Eq. (9). Analytical solution of Eq. (8) gives the exact amplification factor, G r : 


( 10 ) 


G. = c XAt 


Comparing Eqs. (9) and (10), one sees that the three-stage, third-order RK scheme is a polynomial approx- 
imation to the exact solution to third-order. Similarly, the five-stage, fourth-order RK scheme has amplifi- 
cation factor: 


(11) G = 1 + (XAi) + ^(XAt) 2 + ^(XAt) 3 + jjiXAt) 4 + ^(XAf) 

The stability of the RK schemes is shown graphically in Fig. 1 by plotting the |G| = / contour of Eq. 
(9) for the three-stage, third-order scheme and Eq. (11) for the five-stage, fourth-order scheme. A selection 
of XAt which lies in the interior of the closed curve yields \G\ < 1 , he. the scheme is stable. Outside the 
closed curve, |G| > 1 and the scheme is unstable. If the Fourier symbol of the spatial operator, X , is 
purely imaginary (for example the 1-D convection equation) an inspection of Fig. 1 reveals that the region, 
-L f <XAt <Lj , is stable. If X is purely real (for example the 1-D diffusion equation) the region, 
-L n <XAt<0 , is stable. The stability limits for these two extreme cases are given in Table 3 for the 
third- and fourth-order RK schemes. Hence, the fourth-order scheme would allow time steps roughly twice 
that of the third-order scheme. 

Temporal stability analysis of the Navier-Stokes equations is modeled after the stability of the convec- 
tion-diffusion equation for which X is, in general, complex. 

Table 3. Stability limits of Runge-Kutta schemes for purely imaginary (Lj) or real (L R ) spatial operators 


Spatial Operator 

third-order, three stage 

fourth-order, five stage 

Imag, L[ 

1.73 

3.34 

Real, Lp 

2.51 

4.65 


2.3. Spatial Discretization. The numerical approximations to the spatial derivatives appearing in 
the momentum equations, Eq. (3), are given in this section. Standard second-order finite difference approx- 
imations to first derivative terms suffer from large dispersion errors. Spectral methods offer exact differen- 
tiation for resolved modes but suffer from high cost and low flexibility in that simple domains and 
boundary conditions are required for their implementation. In this study, high-order compact finite differ- 
ences are preferred duo to the combination of high-accuracy, flexibility, and relative operation count. 
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2.3.1. First Derivative Terms. Tile first derivative terms appearing in the governing equations are 
approximated using fourth- and sixth-order compact finite difference schemes described by Lele [1], Higher 
accuracy derives from the implicit treatment of derivatives, via: 


(12) 


« 1 >', -I + <!>', + «<!>', + / = + + + 


(or in matrix form: AJ' = B x $ or f = aJ B x + ), where Ax = L X /(N X -1) , N x is the number of 
grid points, <|) ■ represents the first derivative of the generic variable with respect to x , and <x , a, b are 
the coefficients of the compact scheme which determine the accuracy. Similar expressions are used for 
derivatives with respect to the y and z directions. For the fourth-order scheme: a = 1/4 a - 3/2 
and b = 0 , and for the sixth-order scheme: a = 1/3 , a = 14/9 , and b = 1/9 . The LHS of Eq. 
(12) contains the unknown derivatives at grid points i and i ± 1 while the RHS contains the known func- 
tional values <t> ? at the grid points i ± 1 and i ± 2 . A x is a tridiagonal N x x N x matrix and B x is a trid- 
iagonal matrix for the fourth-order scheme and pentadiagonal matrix for the sixth-order scheme. In 
contrast, A x will be diagonal in an explicit finite-difference scheme. 

A comparison of explicit central difference and implicit compact approximations to the first derivative 
is given in Table 4. It is seen that the implicit treatment of the derivative allows for more "compact” sten- 
cils, for given order. Also, the coefficient of the leading truncation-error term is reduced by a factor of 4 for 
the fourth-order scheme and 9 for the sixth-order scheme in comparison to explicit central difference 
schemes of the same order. 


Table 4. Comparison of explicit central difference and implicit compact approximations to the first 

derivative 


Scheme 

Truncation error 

Stencil Size 

fourth-order central 

(-4/5!)( Ax)V 5) 

5 

fourth-order compact 

(-//5!)(A:e)V 5) 

3 

sixth-order central 

(-36/7!)(Ax)V 7) 

7 

sixth-order compact 

(-4/7!)(Aa:)V 7) 

5 


The set of Eq. (13) at all grid points results in a tridiagonal system of algebraic equations and that is 
solved efficiently by factoring the LHS into a lower/upper (LU) system once at the beginning of the simu- 
lation. The LU factors are stored and then used to solve Eq. (12) for the unknown derivatives. 

The resolution properties of the numerical approximation to the first derivative can be analyzed by 
transforming the 1-D convection equation from physical to wavenumber space[l]. For explicit finite-differ- 
ence schemes, a = 0, the numerical wavenumber is given by: 


(13) 


N 


~ Ax X a i' 


i l kAx 


l = -N 
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while 1 for the tridiagonal compact scheme, the numerical wavenumber is given by: 


a s in ( k Ax ) + - cos ( 2 k Ax ) 

1 + 2acos(kAx) 

Note that in general the numerical wavenumber, k*, is complex while the exact wavenumber, k, is real. 
For the 1 numerical approximation to yield an exact solution, the following two conditions must be met: 

(15) Real(k ) = k 

* 

(16) Imag(k ) = 0 

Deviations from Eq. (15) indicate dispersion errors due to odd derivative terms appearing in the truncation 
error, and deviations from Eq. (16) indicate dissipation errors due to even derivative terms appearing in the 
truncation error. The real and imaginary parts of Eqs. (13) and (14) are plotted separately in Fig. 2, for 
several differencing schemes. All four approximations do a reasonable job of approximating the exact wave- 
number (i.e. very low dispersion errors) at, very low wavenumbers ( kAx —> 0 ), but they all do a poor job 
at very high wavenumbers ( kAx —> n ). For intermediate wavenumbers, the fourth- and sixth-order com- 
pact schemes provide a. much better approximation to the exact wavenumber over a greater range of wave- 
numbers than the explicit schemes. The second-order central difference scheme yields a poor approximation 
to the exact wavenumber for all but the very lowest wavenumbers ( kAx <0.5 ). Due to symmetry, central 
difference schemes always have real wavenumbers, hence they contain no dissipation errors. Only the third- 
order upwind scheme has numerical dissipation errors. Spectral methods yield exact differentiation for all 
modes which can be resolved on the specified grid and thus correspond to the exact relationship for kAx in 

Fig- 2. 

Table 5 lists some quantitative measures of resolution for the five schemes. The wavenumber, k ( . y 
defines the region of acceptable accuracy, i.e. for 0 < k < , \k Ax — /rAr| < 0.01 . Modes with k > k c , are 

not accurately resolved. The quantity, k * nax , defines the maximum value for the modified wavenumber, i.e. 
for k<k* mar , the slope of the curve is zero. Also listed is the number of spatial grids points per wave- 
length, PPW = 2n/(k*Ax) , to accurately resolve a given mode. From the estimate of PPW, roughly 
five times as many points are required for the second-order central difference scheme to achieve the same 
accuracy as the compact schemes. 
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Table 5. Resolution measures of various numerical approximations to the first derivative 


Spatial Scheme 

k c Ax 

k maA x 

Points per 
wavelength 

second-order central 

0.22 

1.00 

28.6 

third-order upwind 

0.44 

1.27 

14.3 

fourth-order compact 

1.11 

1.73 

5.6 

sixth-order compact 

1.55 

2.00 

4.1 

spectral 

71 

71 

2 


For 11011 -periodic boundaries, one-sided finite difference expressions are required to close the system of 
equations at the boundary points; i = 1 and i = TV for the fourth-order scheme and i = 1, 2 and i = TV- 1, 
A for the sixth-order scheme. A third-order compact boundary scheme is used at i — 1 and i = TV with the 
fourth-order interior scheme: 


(17) 


*'/ + “6,* 2 = £ X 


I = 1 


where a bs — 2 and 5/2, a bs ^ = 2, - 1/2 are the coefficients of the third-order boundary 

scheme. A similar equation is used at i = TV. 

For the sixth-order scheme, a boundary and near-boundary scheme are required for closure since the 
interior stencil is pentad iagonal. A fifth-order explicit boundary scheme is used at points, i =1 and i = TV: 


(18) 


with coefficients: 


O'; 

©7 

-o 

H > 

II 

a, = -296/105 

a bs s = -215/12 

a bs 2 = 415/48 

a bs 6 = 791/80 

a. = -125/8 

0b 3 

a. = -25/8 

os 7 

a b, 4 = 9 ^/48 

a bs H = 1457336 

near boundary scheme 

is used at points, i 
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( 19 ) *2 = 

i = / 

with coefficients: 

a nb l = a„ 6j = 115/144 

% l>2 = -2U/180 a nl>6 = -7/5 

= 109/48 “nl, 7 = 23/240 

<>„,, 4 = ~ 35/24 = - //72 

The boundary schemes given by Eqs (18) and (19) were shown to be asymptotically stable by Carpenter 
[15]. Similar equations for the boundary and near boundary schemes are used at points i — A and i — N-l. 

2.3.2. Second Derivative Terms. The second derivative terms present in the viscous terms of the 
moment um equation and the Laplacian operator of the Poisson equation for pressure are approximated 
using fourth- and sixth-order compact finite differences. Again, higher accuracy is achieved by treating the 
derivative 1 implicitly: 


a<r, _ / + ♦' + «<t>", + / = - JL - 2 ^i + r 2 4 , i + ♦.■- 7 ) 

(Ax) 

(20) +— ^l(<t>,; + 2- 2< t>, + < t b-2> 

4( Ax) 

(or in matrix form: A xx ty" = B xx ( t> or fy" = A xx B xx ty ), where (b t represents the second derivative of 
the generic variable 4> ( with respect to x, and o i, a, b are the coefficients of the compact scheme. For the 
fourth-order scheme: a = 1/10 , a = 6/5 , and b = 0 , and for the sixth-order scheme: a = 2/11 , 

a = 12/11 , and b = 3/11 . The tridiagonal system of algebraic equations for the second derivatives is 
solved for in the same manner as the first derivatives. A comparison of explicit central difference and 
implicit compact approximations to the second derivative is given in Table 6. As with the first derivative, 
the implicit treatment of the second derivative results in a smaller stencil size for a given order. The lead- 
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ing truncation error term for the compact formulation is reduced by factors of 2 for the fourth-orde scheme 
and 4 for the sixth-order scheme, in comparison to explicit central difference schemes of the same order. 

Table G. Comparison of explicit central difference and implicit compact approximations of the second 

derivative 


Scheme 

Truncation error 

Stencil Size 

fourth-order central 

(~8/6!)( AaoV 6) 

5 

fourth-order compact 

(-3.6/6 ! )( Ax) V 6> 

3 

sixth-order central 

(-72/5! )( Ax)V S) 

7 

sixth-order compact 

(~16.7/8\)( Ax)V S) 

5 


For non-periodic boundaries, one-sided finite differences are required to close the system of equations. 
At i = 1 and i = N, a third-order compact boundary scheme is used: 


( 21 ) 


l = 1 


where a bs — 11 and a bs * — 13, a bs ^ 27, a bs * — 15, and a b = -1 are the coefficients of the third- 

order boundary scheme. For the sixth-order scheme, a near boundary scheme is required at i = 2 and i = 
N-l. The fourth-order interior scheme is used at these points since only a three- point stencil is needed. 

A similar analysis of the 1-D diffusion equation is used to investigate the resolution qualities of the pro- 
posed compact approximation to the second derivative. The tridiagonal compact scheme represented by 
Eq. (20) yields a numerical wavenumber: 


( 22 ) 


♦ 2 
(k ) = 


(AxY 


2a[l - cos(kA: r)] + -[/ - cos{2k,Ax)] 
1 + 2ctcos(kAx) 


The dissipation errors for the explicit second-order central difference, and fourth- and sixth-order com- 
pact schemes are shown in Fig. 3. We see that the numerical wavenumber for the compact scheme more 
closely approximates the exact wavenumber over a wider range of wavenumbers. Quantitative measures of 
resolution power for the various schemes are given in Table 7. Estimates of the PPW show that roughly 
twice as many points are required for the explicit second-order central difference to produce the same accu- 
racy as the compact scheme. 
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Table 7. Resolution measures of various numerical approximations to the second derivative 


Spatial Scheme 

e*A.r) 

( k max Ax ) 2 

Points per 
wavelength 

second -order central 

0.57 

4.00 

11.0 

fourth-order compact 

1.14 

0.00 

5.5 

sixth-order compact 

1.52 

6.80 

4.1 

spectral 

n 

2 

K 

2 


2.4. Enforcement of the Continuity Equation and Poisson Equation for Pressure. An 

examination of the governing equations reveals four scalar equations (continuity and three scalar compo- 
nents of the momentum equation) in terms of four unknowns (three velocity components and pressure). 
Time derivatives for the velocity components in the momentum expiation are used to march those expia- 
tions in time. However, in an incompressible flow, no such time derivative exists for pressure, but the con- 
tinuity expiation gives an additional constraint on the velocity field, in order to enforce the divergence- free 
condition. The current approach overcomes this problem by taking the numerical divergence of the dis- 
cretized momentum equation and substituting in the discretized continuity equation. This results in a Pois- 
son equation for pressure which is solved to ensure that the velocity field is divergence free at the M+l th 
substage. The problem does not arise' in a compressible flow since the density appears as a natural choice 
for the fourth variable, and the continuity equation contains its time derivative to be used for time- 
advancement. The pressure is then obtained from an equation of state. 

Applying the divergence operator S r to the discretized momentum Eq. (3) gives: 


(23) 


b M+1 At 


8. r (^ ) 


= 8 


//; U -8 r 8 x F u 


The term 8 u + , represents the discretized continuity equation at the M 4- 1 sub-stage and is set to 
5 

zero to enforce the divergence-free condition. The term, 8 X u\ , represents the divergence of the velocity 
field at the previous substage M. This term should also be zero, but in practice, it is retained to “kill off’ 
any accumulation from previous substages due to lack of convergence, etc. The term, 8 x H] , is the source 
term of the Poisson equation and represents gradients of the convection and diffusion terms which are 
known from the previous sulvstage. The term, 8 X 8 x P" , represents the discretized Laplacian operator on 
the pressure. Hence, the Poisson equation to be solved for the pressure is: 


(24) 


vV W = 8,. 


M 


H M + 


b Al+I At 


The solution details follow. 
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3. Solution of the Poisson equation for pressure. A significant amount of the total computa- 
tional time required for the solution of the incompressible Navier-Stokes equations (as much as half) is 
devoted to the enforcement of the continuity equation/solution of pressure. This stems from the fact that 
evolution equations exist for the velocity components (i.e. the momentum equations) while none exist for 
the pressure. Instead, an elliptic equation must be solved for the pressure which involves the solution of a 
system of equations and is expensive. Solutions methods for elliptic equations generally fall into two catego- 
ries - direct or iterative. Direct methods usually involve some form of Gaussian elimination where the coef- 
ficient matrix is first factored into an upper and lower matrix and then the solution is computed using back 
substitution. Tin? operation count and memory requirements for this procedure can be prohibitively large 
for the solution of systems involving a large number of unknowns ( - JO 6 in typical 3-D problems). The 

alternative to a direct solution is an iterative procedure where an initial approximation to the solution is 
used to yield an improved solution. This process is repeated until the solution is converged to within a 
desired tolerance. The operation count and memory requirements of most iterative methods are less than 
that of Gaussian elimination. Therefore, the iterative solution procedure is used in this study to solve the 
Poisson equation for pressure. The details of this procedure are outlined in this chapter. 

3'1. Discretized Laplacian Operator. The discrete Poisson equation for pressure which was 
derived in Sect. 2.4 is given by: 


(25) 


M, '' 


M 


= 5 


M 


h m + 


,M + 1 ., 
b At 


The LHS of Eq. (25) represents a discretized Laplacian operator composed of two applications of the 
first derivative operator, 8 r . It. is well known that using two first derivative operators to represent the 
Laplacian operator on non-staggered grids can lead to an “odd-even” decoupling of the solution. Indeed, 
with standard second-order central differencing for the first derivative operator, the solution at even grid 
points completely decouples from the odd grid points, leading to non-physical results. One remedy is to 
introduce terms of the same order as the truncation error which in effect replaces the two first derivative 
operators with a single second derivative operator. This couples the solution at odd and even grid points 
while maintaining the same formal order of accuracy. The Laplacian operator is discretized using a single 
second derivative operator to prevent possible decoupling and Eq. (25) becomes: 


(26) 


8 . = 8 


M 


h m + 


b M + l At 


where S XX P^ represents the discrete Laplacian of pressure and is discretized using the compact second 
derivative operator given by Eq. (20). 

The effect of replacing the two first derivative operators with the single second derivative Laplacian 
operator is shown by rewriting Eq. (26) as 
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M 


(27) 


8,,r r U = 8., 




b A,+ l At 


-R 


where an iterative solution is assumed and P is the current iterate of pressure and R is the residual 
imbalance due to incomplete convergence. Taking the divergence of the momentum equation Eq. (3) after 
the iterative solution of the pressure and time advancement gives: 


M 


(28) 


5 2 p u = 5 , 


h A! + 


b M + , At 


- 8 ,. 


M+ 1 



M + 1 

yu OIJ 


b At) 

Subtracting Eqs. (27) from (28), and rearranging gives the expression for the divergence: 


( 2 !)) 


8 ,. u 


A1+ 1 


= -b AI + , At\K -8 


p u + r 


Equation (29) shows that the imbalance of the continuity equation is due to two contributions; (l) the 
difference lietween the two different Laplacian operators which is ()(A.r ) where m is greater than the 
order of accuracy of the Laplacian operators (it is easy to show that it is fourth-order for explicit second- 
order operators) , and (ii) the residual imbalance due to incomplete iterative convergence. Equation (29) 
also shows that the residual imbalance, R, does not need to be driven to machine zero, only below the level 
of the difference in the two Laplacian operators. However, the convergence tolerance for terminating the 
iterative solution of the Poisson equation must be sufficiently small for the solution to be insensitive to the 
stopping criterion. As the grid is refined, the tolerance must be reduced in accordance with the order of 
accuracy of the scheme. For example, when using the sixth-order scheme on relatively fine grids (see Table 
18), it is necessary to converge the solution to machine zero (1()" 14 ) in order for the solution errors to reduce 
at a sixth-order rate. 

Equation (26) at all interior grid points results in a system of equations that is solved at each sub-stage 
of the time advancement scheme. For simplicity, the system of equations is described for the 2-D Poisson 
equation with periodic boundaries on a uniform grid, for the fourth-order tridiagonal scheme defined in Sec. 
2.3. The solution procedure is easily extended to the 3-D Poisson equation, as presented in Api>endix A. 
Extension to general curvilinear grids is presented in Appendix B. 

The set of Eq. (2G) can be written in the form of a system of equations as: 


(30) 


ap= [A Lb 
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(nj x nj block matrix) (ni x ni scalar matrix) 


I a7 a7 

al 1 a I 

o 

Ayy = 0 , where I is the ni x ni identity matrix 

o 

a! I al 
Oi l a l I 

(nj x nj block matrix) 



(nj x nj block matrix) (ni x ni scalar matrix) 



(nj x nj block matrix) 


P= [ P/ p 2 . o o pj r . 

where - 

[ P U P 2,J ° ° ° p n J 

(nj x 1 block vector) 

(ni x 1 scalar vector) 

F = \Fi F 2 °°° F n f, 

where Fj = j 

[ F U F 2,J ° ° ° F m,l 

(nj x 1 block vector) 

(ni x 1 scalar vector) 


where P ui and F ifj = 5 x [H i + u/(b‘ ! + 1 M)]^ . are the pressure and source terra at the i,j grid point, 
respectively. The symbols, m = N x +1, nj = N y +1 denote the number of grid points in the x, y directions, 
respectively. Values at i = N x +1 and j = A^+l locations are replaced by values at i = 1 and j = 1 since 
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the boundaries are periodic. This change results in a lion-zero element in the upper right and lower left cor- 
ners of the coefficient matrices. The sixth-order compact scheme (with pentadiagonal RHS) can be written 
in the same manner by including the i±2 and j±2 terms in tin' B XY and Byy matrices. 

For non-periodic boundaries, the second derivative boundary scheme given by Eq. (21) is used at the 
txnmdary points. In addition, tin? unknown pressures at the boundaries are replaced with their boundary 
values and those terms are moved and added to the RHS of Eq. (26). For Dirichlet boundary conditions, 
such as freest ream conditions, this procedure is straightforward. Neumann boundary conditions, such as 
those applied at inflow and outflow planes, require that the pressure gradient at the boundary be dis- 
cretized using a first derivative scheme (Eq. 17): 

3 

(si) = 

i = / 

when- the subscripts “l,j” art- used to denote the inflow plane for example. The boundary pressure, Pjj is 
then solved for: 


(32) 


3 


p u = 



'L a b» i p i,j 

i = 2 


Equation (32) is then used to substitute for the boundary pressures in Eq. (30). Tin- first term on the RHS 
of Eq. (32) is known from the boundary condition and is moved and added to the RHS of Eq. (30). The 
second term on the RHS of Eq. (32) contains the unknown pressures, P 2 j and P 3j , so they are kept on the 
LHS of Eq. (30) and corresponding terms of the coefficient matrix, A are modified accordingly. The result- 
ing system of equations contains only the interior pressures as unknowns, P, r 2 < i < N x - 1 and 

2<j<N, r l . 

Equation (30) results in a “cross” type stencil at the i, j node in which all points along lines passing 
through the central node contribute to the stencil. The coefficients of this stencil are implicitly defined in 
the sense that the matrix operations, {A~J g B xx + A~ y l y B yy ] , in Eq. (30) must he performed to determine 
their values. 

Multiplying Eq. (30) by A yy A xx gives: 


(33) \A yy D xx + A yy A xx A' y ‘ y B yy )P = A yy A xx F 

It is easy to show by inspection that matrices A xx and A yv commute, i- e - AyyA xx A xx A yy , so that Eq. 

(33) simplifies to: 

(34) l A yy B xa- + A xr B yy\ P = A yy A x. r F 

With the coefficients of the fourth-order compact approximation, the second derivative results in an 
explicit nine-point, “grid” type stencil for the LHS and R.HS of Eq. (34). 
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With the coefficients of the sixth-order compact approximation, tile second derivative results in the 
same: nine-point stencil for the RHS and an explicit twenty-one point stencil on the LHS: 


(36) 
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( P I +2,j + ' P i -2,j] + 

[ P i+2,j+l + Pi + 2,j-l + Pi -2, j + 1 + Pi-2, j-l] + 

I P i + l,j + 2 + P i + I,j-2 + P i-I,j + 2 + Pj-J,j-2^ + 

l P i + I,j + 1 + P i + l,j- 1 + Pi- l,j+ 1 + Pj- l,j - il = F i,j + 


a[F iJ + i + F i.j-l + F i + l,j + F i-l, j ) + <* 2 lF l + ,,j +l + F i + l j - 1 + F l _ ]j + 1 + F l _ Ij _ I ] 


Thus for uniform cartesian grids, the stencils of Eqs. (35) and (36) are preferred because they require 
fewer operations and their coefficients are explicit in contrast to the stencil of Eq. (30). The commutative 
property of the A xx and A yy matrices is valid even with non-periodic boundaries. Numerical experiments 
confirm that the cross-type stencil represented by Eq. (30) and the grid-type stencil of Eq. (34) give identi- 
cal results. 
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3.2. Iteration Matrix. The system of equations represented by Eq. (30) or Eq. (34) must be solved 
to convergence to ensure a divergence- free velocity field, at every sub-stage of every time step. Although 
the matrix is sparse, direct sparse-mat rix methods were' found to be competitive only for fairly coarse-grid 
2-D problems Therefore, iterative schemes are required, and point-relaxation methods are utilized. In this 
approach, only the value at the central node of the stencil, P iJt is treated as an unknown so that the multi- 
diagonal system of equations degenerates to a diagonal system for one relaxation sweep, which is trivial to 
solve. This process can be written in matrix notation bv decomposing the matrix, A, into the sum of the 
diagonal, lower, and upper matrices of A: 


( 37 : 


AP = [D- L — U]P = F 


where the matrix, D, is the diagonal matrix of A, and the matrices, L, U, are the lower and upper matrices 
of A, respectively. The solution at the current iteration level, P*, is corrected with the increment, P' , to 
yield the solution at the next iteration level, P = P* + P f - For weighted Jacobi iteration the pressure is 
calculated from: 


(38) P = P + wD P 

where 0< o> < 1 denotes under-relaxation. Thus, Jacobi iteration is equivalent to computing the residual 
(ft* _ p _ a P*) of the current iterate, P*, at all grid points followed by an update operation. In this 
regard, information is held and the solution is updated at all grid points simultaneously. Since the compu- 
tation of the residual vector and the update of the solution vector are completely separate operations, eac h 
operation is fully vectorizable. This results in an improved computational rate' when those operations art' 
performed on vector computers. 

On the other-hand, weighted Gauss-Seidel calculates the pressure from: 


(39) P = P* + w[D-L] ‘ R 

which requires that the solution vector be updated in ascending order (P\.\, -Pj.n— P ni,nj )• Since this does 
not vectorize very well red-black ordering is a possibility. Although the Jacobi method converged more 
slowly than the Gauss-Seidel method, the computational speed increase due to vectorization led us to pre- 
fer it in the present study. 


3.3. Multigrid Solution. Relaxation schemes such as Jacobi and Gauss-Seidel applied on a single 
grid level suffer from poor convergence rates as the number of grid points increase. Multigrid methods over- 
come these deficiencies by utilizing a hierarchy of grids. Smooth error components are transferred to 
coarser grids where they appear as high frequency error components and are quickly removed by relaxation 
sweeps. A coarse grid correction scheme is utilized in the current study to improve the convergence rate of 
the pointwise relaxation scheme on a single grid. Subscripts are used to denote grid level, i.e. Ph and P ah 
denotes the solution on the fine and coarse grids, respectively. The symbol, I h , is used to denote transfer 
from the fine to the coarse grid, while j£ ; . is used to denote transfer in the opposite direction. The algo- 
rithm for one coarse grid correction is given below and additional details can be found in Briggs [16]. 

(1) Smooth the current iterate, P* h , on the fine grid times: 
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A h^ = F h 


(2) Calculate the residual on the fine grid: 


(3) Transfer (restrict) tile residual to tile coarse grid, where it is used as the source train for the error 
equation: 


F - J 2il /? 

r 2h “ J h K h 


(4) Solve for the error on the coarse grid: 


E 2„ = A 2k F 2k 


(j) Transfer (prolongate) the error to the fine grid and correct the solution: 


v 


(6) Perform v 2 post-relaxation sweeps: 


= F „ 


Standard second-order interpolation is used to transfer variables from the coarse to the fine grid, while 
the full weighting operator [10] is used to transfer variables in the opposite direction. Although the above 
algorithm utilizes only two grid levels, improved efficiency results from incorporating as many grid levels as 
possible. In this respect, the direct solution of the error equation in step (5) is performed on a very coarse 
grid requiring a small number of operations. 

Since the simulations are performed on vector computers, Jacobi iteration was utilized for relaxation 


sweeps because it is fully vectorizable. Two pre- and two post-relaxation sweeps were performed on each 

grid. Through numerical experiments, the optimum relaxation factor for the uniform grid formulation was 
found to be, co = 0.9 


Convergence rates (CR) for the 2D wave decay problem (discussed in Sect. 4) using the 4th- and 6th- 
order schemes are given in Table 8 along with the number of multigrid levels used. The results demonstrate 
that the convergence rate of the multigrid solution technique is independent of grid spacing unlike typical 
single grid, iterative solution techniques where the convergence rate approaches unity as the grid is refined. 

Effects of grid aspect ratios on multigrid convergence rates are shown in Table 9 for the solution of the 
2D Stuart’s problem using the fourth-order compact scheme. The 2D Stuarts problem, which is defined in 
Sect. 4, is a temporally developing mixing layer with neutral growth, which is governed by the 2D Euler 
equations, and for which an exact solution exists. Uniform grid spacing is used in the a: coordinate direc- 
tion, while a logarithmic function is used to cluster grid points in the y coordinate direction where large 
gradients exist due to the mixing of the high and low speed freestreams (see Figure 10): 
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(40) 



where r is a clustering parameter and T) is the uniformly distributed computational coordinate. The results 
in Table 9 show that the multigrid convergence rate deteriorates somewhat with increased grid aspect, 
ratio. 

As more grid points are moved to the mixing region through clustering, solution errors in the y coordi- 
nate direction decrease while those in x direction remain constant. Thus, the overall solution error eventu- 
ally saturates with increased resolution in the y direction 


Table' 8: Multigrid convergence rates, CR, for solution of 2D wave decay. 


ni x nj 

number of MG 

levels 

... 

CR - 4th order 
scheme 

CR - 6th order 
scheme 

16 x 16 

2 

0.41 

0.30 

32 x 32 

3 

0.40 

0.30 

64 x 64 

4 

0.41 

0.30 

128 x 128 

5 

0.42 

0.30 


Table 9. Multigrid convergence rates, CR, for solution of 2D Stuarts problem on non-uniform grids 

(nixnj = 49x 101) at t,=0.2. 


Clustering 
parameter, c 

< < 
II 

CR 

N 2 

uniform grid 

1 

0.19 

0.1319 x 10~ 4 5 

0.12 

3.2 

0.21 

0.7860 x 10 6 

0.16 

6.2 

0.28 

0.6402 x 10‘ 6 

0.20 

12.2 

0.38 

0.5725 x 10" 6 

0.24 

24.2 

0.48 

0.5675 x 10 ° 


4. Results and discussion. The performance of the numerical formulation is tested by application 
to a variety of benchmark problems. Emphasis is placed on the numerical approximation of spatial deriva- 
tives. In particular, the convection terms (containing first, derivatives) present the most difficulty in numer- 
ical approximation since large dispersion errors exist at high wavenumbers ( kAx ~n ). It is essential that 
the numerical approximation to the first derivative provide low dispersion errors over a large range of 
wavenumbers. This is especially true in 3-D simulations where reducing the required number of grid points 

by half in each coordinate direction leads to eight times fewer total grid points, and considerable savings in 

computer time 
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In Section 3, theory showed that compact schemes require roughly five times fewer points to accurately 
resolve a given mode in comparison to the standard second-order central difference approximation to the 
first derivative. This theory is tested by solving some practical problems ranging from the 1-D convection 
equation to the 2-D Navier-Stokes equations. 

4.1. 1-D Convection Equation. The first, problem to be solved is the 1-D convection of a Gaussian 
profile: 


(41) 


du 

di c dx 


= 0 


subject, to: u(x, 0 ) = 0.5 exp -0) In (2) 


-20<x<450 ; Ax = 1 , c = 1. 


This could test the time advancement scheme and the numerical approximation to the first derivative. 


The exact solution corresponds to the convection of the initial profile at the constant wave speed, c. The 
third-order RK scheme was used to advance the equation in time for all spatial schemes. In addition, the 
CFL number was kept small so that resulting errors are due solely to the spatial formulation. Distortion in 
the shape of the profile indicates dissipation and/or dispersion errors in the solution. The convection equa- 
tion was solved using three approximations to the first derivative; (i) a second-order central difference, (ii) 
a third-order upwind, and (iii) the fourth-order compact approximation outlined in Section 2.3. The 


parameters and initial conditions are those proposed at the ICASE/LARC Workshop on Benchmark Prob- 
lems in Computational Aeroacoustics, Hardin et al. [17]. Since the specified grid is relatively coarse, this 
problem provides an excellent test of the resolution power of the numerical approximation. Figure 4 shows 
the computed solutions at t = 400 after the profile has convected to x = 400. There is little discernible dif- 
ference between the exact solution and the solution with the fourth-order compact scheme. However, the 
solutions w r ith the second-order central difference and the third-order upwind approximations show greatly 
reduced peak values and large, dispersive waves trailing the Gaussian profile. The errors from the second- 
order central difference scheme are the most severe. 


It is difficult to determine by inspection what portion of the error is dispersive and w r hat portion is dis- 
sipative. The solutions are transformed into wavenumber space using a Fourier transform method and com- 
pared with the exact solution in Fig. 5 to address this issue. The graph displays the resulting complex 
Fourier coefficient in polar form with the amplitude displayed in Fig. 5a and the phase angle in Fig. 5b. It 
can be seen from Fig. 5a that the solutions computed with the second-order central difference and fourth- 
order compact schemes predict the correct amplitude for all modes. The amplitude of the solution com- 
puted with the third-order upwind scheme is reduced or dissipated, especially at higher wavenumbers. Fig- 
ure 5b shows that the fourth-order compact scheme predicts the correct phase angle even for the highest 
wavenumbers. 

The phase angle from the second- and third-order solutions are only correctly predicted for the very 
lowest wavenumbers ( k<0.2 for the second-order solution and k<03 for the third-order solution). 
Large dispersion errors are evident at high wavenumbers. The above trends in the numerical solutions are 
consistent with the dissipation /dispersion error theory for the 1-D convection equation and show the reso- 
lution pow r er of the compact schemes. 
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The second problem, also proposed at the ICASE/LARC Workshop on Benchmark Problems, is the 
solution of the 1-D convection equation in a spherical coordinate system. In Cartesian coordinates, the gov- 
erning equation takes the form: 

du u du n 

(42) a7 + ; + aI = 0 

subject, to: u(x,0) = 0 ■ u(5,t) = sin(^ — j ; 5 < 450 ; Ax = I 

Figure 6 shows the exact solution at t. = 400 which corresponds to a damped sine wave due to the addition 
of the u/x term in the governing equation. Fig. 7 shows computational results for the region, 200 < x < 220 
using the third-order upwind approximation to the first, derivative on three different grids and the fourth- 
order compact scheme, on the specified grid only, which corresponds to 8 PPW. The upwind solution with 
8 PPW shows severely reduced amplitude and a phase shift, relative to the exact solution. Even those with 
10 PPW and 32 PPW arc' not very accurate. It takes roughly 64 PPW (not shown) t.o reproduce the exact 
solution with the third-order upwind approximation. On the other hand, the fourth-order compact approx- 
imation is able to reproduce the exact solution with 8 PPW . 

4.2. 2-D Convection Equation. Multidimensional effects of the numerical formulation are explored 
by solving for the convection of an inverted cone around a circle. This problem is governed 2-D convection 
equation: 

du 3 u du n 

(43) dt +C ^ +C ^j = ° 

where c x = -y and c = x, are the convection speeds in the x and y directions, respectively. The initial con- 
ditions are that of an inverted sharp cone centered at x, y = -0.5, 0. The exact solution corresponds to the 
cone being converted counterclockwise in a circular path of radius, r 0 = 0.5 with a period of 2n . Distortion 
of the shape of the coni' is an indication of dispersion and/or dissipation errors. 

Figure 8 shows computed results after one revolution of the cone using (a) a third-order upwind 
approximation and (b) a fourth-order compact approximation to the first derivatives on a 32 x 32 grid with 
uniform spacing. This grid defines the shape of the cone with a maximum of 8 points in each coordinate 
direction. The exact shape of the cone is included to the right of the computed solution at x, y = 0.5, 0 for 
comparison purposes. The third-order solution (Fig 8a) shows that the sharp point, of the cone is greatly 
diffused and that dispersion errors are evident trailing the cone. A grid of 128 x 128 (or 32 points defining 
the' shape of the cone) must be used with the third-order upwind approximation before the shape of the 
cone is faithfully reproduced. The fourth-order compact solution (Fig 8b) shows that the shape of the cone 
is not distorted as it is converted around the circle on the 32 x 32 grid. Indeed, the only noticeable error is 
a very small “grid to grid” oscillation due to the absence of physical viscosity in this problem and numerical 
viscosity in the compact scheme. 

Figure 9 shows results for the same problem after one revolution obtained by Orszag [18) using (a) sec- 
ond-order Arakawa finite-difference, (b) fourth-order Arakawa finite-difference, and (c) a spectral approxi- 
mation to the first derivatives on a 32 x 32 grid. The finite difference solutions show errors similar to the 
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third-order solutions in Fig. 8. The spectral method, which provides exact differentiation for all wavenum- 
bers representable on the 32 x 32 grid, convects the cone without distorting its shape. Thus, the solution 
using the compact scheme is closer to the spectral solution than the solutions obtained with conventional 
finite difference schemes. The higher accuracy and resolution characteristics are achieved by the implicit 
treatment of the derivative. Even though the stencil size of the compact scheme is finite, the implicit treat- 
ment of the derivatives makes the scheme global much like that of spectral methods. 

4.3. Euler/Navier-Stokes Equations. In the previous sections, the effect of numerical approxima- 
tion on the accuracy of the convection terms was documented. In this section, the accuracy of the enforce- 
ment of the continuity equation through the solution of the Poisson equation for pressure is documented by 
solving the 2-D Euler/Navier-Stokes equations. Since the Navier-Stokes equations contain viscous terms, 
the numerical approximation to the second derivative is also tested. The test problems chosen for valida- 
tion contain many features of the 3-D jets which are simulated in the current study. In this respect, the 
test problems are not merely academic exercises. There benchmark problems are solved; (i) a temporally- 
developing plane mixing layer (2-D Stuarts problem), and (ii) 2-D viscous wave decay Problems (i) and 
(ii) have exact solutions. 

4.3.1. Temporally-Developing Plane Mixing Layer. Exact solutions to the Euler or Navier- 
Stokes equations for general flows do not exist, due to the non-linearity of the convection terms. However, 
under special conditions exact solutions may be found. An exact solution for the temporally-developing 
mixing layer was first published by Stuart [19], The initial conditions for the 2-D Stuarts problem corre- 
spond to a steady hyperbolic tangent function for the streamwise velocity component with a periodic array 
of vortex cores in the mixing region which cause the solution to vary in time. The wavelength of the distur- 
bance corresponds to the neutral mode such that the disturbance is convected in the streamwise direction 
with no change m amplitude. The exact solution for the streamwise velocity component, u, and the trans- 
verse velocity component, v is given by: 


u(x, y, t) 


c+ Csinh(u) 

Ccosh(y) + Acos(x — ct) 


(44) 


*'(*» y, 0 


Asin(x — ct) 

Ccosh(y) + Acos(x - ct) 


where A = Jc* - I is a parameter which controls the strength of the perturbation and c is the convec- 
tive speed of the mixing layer. The flow is periodic in the streamwise direction with length, L r = 2n , 
0<x<2n .The flow is infinite in the transverse direction but in this study is truncated at a finite dis- 
tance, -L y <y< L y , such that a zero-traction freestream boundary condition is well approximated. Tests 
which vary the transverse domain height, 2 L y , show that L y = 10 is sufficiently large to implement this 
boundary condition. The exact solution is shown in Fig. 10a with parameters, c = 1, A = 1/2. A uniform, 
cartesian grid is used for the simulations in this section. Unless otherwise specified, the third-order RK 

scheme is used for time advancement and time steps are sufficiently small so that spatial errors are dorni- 
nant. 
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Figure 10b shows the numerical solution at. i = 20 k (ten flow through times) on a relatively coarse 
grid of 13 (streamwise) x 41 (transverse) using the fourth-order compact approximation of convection 
terms and pressure. The solution of pressure involves the computation of the source* term and the discreti- 
zation of the Laplacian operator in Eq. (33). In addition, once the Poisson equation is solved for pressure, 
the gradient of pressure is computed which is required to advance the momentum equation in time. There- 
fore, the phrase “fourth-order solution of pressure” corresponds to the source and pressure gradient terms 
computed with the compact first, derivative scheme outlined in Section 2.3.1 while the Laplacian operator is 
discretized using the compact, second derivative scheme* outlined in Section 2.3.2. 

Even though the grid is relatively coarse ( 13 streamwise points per wavelength and roughly 8 points in 
the mixing region at y-0 ), there is little discernible* difference between the exact and numerical solutions 
after ten flow through times. It is important to check the convergence of the* error as the grid is refined to 
expose any coding errors, to demonstrate that the order of error convergence seen in practical computa- 
tions is that predicted by a Taylor series analysis, and to gain confidence in the numerical formulation. 
Tables 10 and 11 give a quantitative measure of the L2 and maximum errors in the velocity components at. 
t — o.i using the fourth- and sixth-order compact, approximations to the convection terms and solution of 
pressure, respectively. Solution errors from three grids are shown where the grid spacing in the x and y 
directions is halved from coarsest to finest, grid. The results show that the L2 and maximum errors con- 
verge at roughly the rate predicted by a Taylor series analysis as the grid is refined. The order, N, is com- 
puted using the solution error from three grids of spacing, h, 2h, and 4h. 


$2h ~ §4h 

In 2 


where + 4k are the errors on the h, 2/i,and 4/i grids, respectively. In using Eq. (45) it is assumed 

that the solution is fully resolved on all three grids and that the leading truncation error term is dominant 
(Demuren and Wilson [20]). 


Table 10. Solution errors for 2-D Stuarts Problem at t = 0.1 using fourth-order compact approximation for 

convection terms and solution of pressure. 



To address the effect of computing the pressure with a lower-order formulation, the 2-D Stuarts prolv 
lem was solved using second-order central, fourth-order compact and sixth-order compact approximation of 
the convection terms but a second-order central difference solution of the pressure. The results of the three 
computations are shown in Table 12. These results show that the lower-order solution of pressure results in 
the overall convergence of the error being second-order, even if the convection terms receive a higher-order 
treatment. All terms must be discretized using higher-order approximations to achieve higher-order error 
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Table 11. Solution errors for 2-D Stuarts Problem at t = 0.1 using sixth-order compact approximation for 

convection terms and solution of pressure. 


Grid(ni x nj) 

Max U Error 

Max V Error 

L2 Norm U 

L2 Norm V 

13 x 41 

0.73 x 1R 3 

o.io x itr 2 

0.97 x 10 4 

0.11 X 10' 3 

25 x 81 

0.17 x 1R 4 

0.20 X 1R 4 

0.12 x 10’ 5 

0.14 x 10 5 

49 x 161 

0.15 x 10 -5 

0.18 x 10 -5 

0.25 x 1R 6 

0.45 x 10' 6 

Order, N 

5.4 

5.7 

6.3 

6.2 


convergence rates. Thus, formulations presented in the literature, such as by Najjar and Taft.i [21], which 
use higher-order differences for the convection and diffusion terms but second-order differences for the Pois- 
son equation for pressure would only be globally second-order accurate. 


Table 12. Errors with different discretization approximations for convection, but second-order discretization 

for pressure. 



Second-order central 

Fourth-order compact 

Sixth-order compact 

Grid (ni x ni) 

L2 Norm U 

L2 Norm V 

L2 Norm U 

L2 Norm V 

L2 Norm U 

L2 Norm V 

13 x 41 

0.21 x 10' 2 

0.20 x 1R 2 

0.20 x 10 -2 

0.20 x 1R 2 

0.20 x 10' 2 

0.20 x 10 -2 

25 x 81 

0.45 x 10‘ 3 

0.53 x 1R 3 

0.42 x 10‘ 3 

0.44 x 10' 3 

0.42 x 1R 3 

0.44 x lO' 3 

49 x 161 

0.11 x 10' 3 

0.13 x 10 -3 

0.11 x 1R 3 

0.11 x 1R 3 

0.10 x 10‘ 3 

0.11 x 10‘ 3 

Order (N) 

2.3 

2.1 

2.4 

2.2 

2.3 

2.2 


.The solution of the 2-D Stuarts problem validates the numerical formulation for the enforcement of 
the continuity equation and the solution of the Poisson equation for pressure. In addition, it has been 
shown that the zero-traction freestream boundary condition for shear flows is a valid approximation pro- 
vided that the freestream boundary is located a sufficiently far distance from the mixing region. 

Temporal accuracy of the overall RK schemes was confirmed by performing computations for the 2-D 
Stuarts problem on the 49 by 161 grid and an even finer 97 by 161 grid, with several different time steps. 
The largest time step in each case was chosen from stability considerations. The results are shown in 
Tables 13 - 15. In Tables 13 and 14, the velocity field is specified from the exact solution and the vorticity 
is computed as a passive scalar, using the sixth-order compact scheme. In each case, the error shown con- 
tains both temporal and spatial discretization errors. As the time step is reduced to very small values the 
latter become dominant and error reduction stagnate as seen in column 5 of Table 13. On this grid, the 
spatial discretization error can be estimated to be about 0.6 x 10' G . On the 97 by 161 grid, the spatial dis- 
cretization error should reduce by about 2 G or 64, to roughly 0.9 x 1R 9 . Column 5 of Table 14 does show 
stagnation at about this value, which confirms the analysis. Equation (45) enables separation of temporal 
and spatial discretization errors by comparison of computed results from three time steps on the same grid, 
or three grids with the same time step. This procedure shows (column 4) that the 3-3 RK scheme is indeed 
third-order accurate, and the 5-4 RK scheme (column 7) is somewhat better than fourth-order accurate. 

The use of computed results from only two time steps (columns 3 and 6) is clearly erroneous and produces 
unreliable conclusions. 
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Table 13. Solution errors for 2-D Stuarts Problem at t, = 1.0, with different RK schemes on a 49 by 161 

grid. 



3-5 RK scheme 

5-4 RK scheme 


Time 

Step 

L2 Norm of 
Vortieity Error 

N from 2 
grids 

N from 3 
grids 

L2 Norm of 
Vortieity Error 

N from 2 
grids 

N from 3 
grids 

0.10 




0.32980 x 10 -5 



0.05 

0.12803 x 10‘ 4 



0.69910 x Ht 6 

2.24 


0.025 

0.17157 x 10‘ 5 

2.9 


0.63575 x 10 ° 

0.14 

5.4 

0.0125 

0.66280 x 10 _c 

1.4 

3.4 

0.63328 x 10 ° 

0.01 

4.7 


Table 14. Solution errors for 2-D Stuarts Problem at t = 1.0, with different RK schemes on a 9i by 161 

grid. 



3-3 RK scheme 

5-4 RK scheme 


Time 

L2 Norm of 

N from 2 

N from 3 

L2 Norm of 

N from 2 

N from 3 

Step 

Vortieity Error 

grids 

grids 

Vortieity Error 

grids 

grids 

0.05 




0.1965 x 10 -6 



0.025 

0.1613 x 10' 5 



0.1683 x 10' 7 

3.5 


0.0125 

0.2019 x 10' 6 

3.0 


0.9496 x 10~ 8 

0.83 

4.6 

0.00625 

0.2685 x 10" 7 

2.9 

3.0 

0.9310 x 10‘ 8 

0.02 

5.3 


In Table 15, the full set of Euler equations is solved with the fourth-order compact scheme and tin' 3-3 
RK scheme. Third-order temporal accuracy is confirmed. In order to achieve this, it, was necessary to 
obtain a divergence-free velocity field at every sub-stage of the time-stepping proc ess. 


Tabic' 15. Solution errors for 2-D Stuarts Problem at t. = 1.0 using third-order RK and fourth-order 

compact approximation for Euler Equations. 


Grid (ni x nj) 

Time Step 

L2 Norm U 

L2 Norm V 

49 x 161 

0.05 

0.51630 x 10 3 

0.25821 x 10~ 3 

49 x 161 

0.025 

0.27651 x 10' 5 

0.18778 x 10' 4 

49 x 161 

0.0125 

0.24666x 10‘ 6 

0.17823 x 10' 6 

Order, N 


3.0 

2.9 


4.3.2. Viscous Wave Decay. The numerical treatment of viscous terms is validated by solving the 
2-D viscous wave decay problem which is governed bv the Navier-Stokes equations. The domain for this 
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problem is periodic in both the x and y directions where periodic boundary conditions are applied. The 
exact solution is given by: 


u(x, y,t) = -cos(x)sin(y)e 


v(x,y,t) = sin(x)cos(y)c ' Re ' 

where Re = 20, L x = L y = 1, and t = 0 gives the initial conditions. The exact solution consists of sinusoi- 
dal waves in the x and y directions which decay in time. Table 16 shows the L2 norm of the error at t = 
0.025 using the fourth- and sixth-order compact approximations for convection and diffusion terms and the 
solution of pressure. The results are compared to the fourth-order, Essentially Non-Oscillatory (ENO) 
scheme from Weinan and Slm[22]. The error converges at fourth- and sixth-order rates thus validating the 
numerical treatment of the viscous terms and again validating the convection terms and the solution of 
pressure. The error of the ENO scheme converges at a fourth-order rate, but is more than two orders of 
magnitude greater than the fourth-order compact results. The error magnitude of the sixth-order compact 
formulation on the 128 x 128 grid has reached the round-off error level ( - ] 0~ 13 ) of the Cray supercom- 

puter, indicating that extremely accurate results are obtained on average-sized grids. 


Table 16. Solution errors for 2-D viscous wave decay. 


Grid (ni x nj) 

4th oa compact 

6th oa compact 

3rd (4th) oa ENO 

16 x 16 

0.14 x 10' 1 2 3 4 * 6 

0.10 X 10' 7 

_ 

32 x 32 

0.77 x 10' 8 

0.15 x 10' 9 

0.53 x 10~ 5 

64 x 64 

0.47 x 10" 9 

0.27 x HT 11 

0.32 x 10~ 6 

128 x 128 

0.71 x 10' 10 

0.11 x 10 12 

0.20 x 10 -7 

Order, N 

4.0 

6.0 

4.0 
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APPENDIX A - Discretization of 3D Poisson equation. In Sect. 3, the solution of the 2D 
Poisson equation for pressure using fourth-order compact finite differencing was developed and a reduced, 
nine-point, grid-type stencil was presented. In this appendix, the details of the solution procedure for the 
3D Poisson equation on uniform grids are presented. The Poisson equation for pressure was given by Eq. 
(26): 


X nM X rjA/ U t 

0 ry F = 5, H- +-*77 — ; — 

xx i x , i ,M+/ 4l 

0 At 


Compact finite differencing is used for the discrete derivative operators in Eq. (A.l). For the 3D Pois- 
son equation on uniform grids, Eq. (A.l) can be written in the form of a system of equations as: 


( A.2) AP = [aJ t B xx + AyyB yy + A z [b zz ]P = F 

where the tilde is used to denote 3D block matrices which are defined using the 2D block matrices given in 
Sect. 3. 
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-21 


-21 1 
I -21 I 


I -21 I 
-21 I -21 

where the above are nk x nk block matrices, nk is equal to the number of grid points in the z direction, and 
I is the (ni x nj ) x (ni x nj) identity matrix. The 2D block matrices A xx , A yy , B xx , and B yy were defined 
in Sect. 3. The block vectors, P and F are assembled from the 2D vectors, P and F, defined in Sect. 3, 
p - J p f P 2 ° ° ° F = | F } F 2 ° ° ° F, ?A J . The sixth-order compact scheme (with pentadiagonal RHS) 

c*an be written in the same manner by including the i±2 , j±2 , and k±2 terms in the above matrices. The 
coefficients of this stencil are implicitly defined in the sense that the matrix operation, 

[A~j! x B rr + Ay\ ( B uu + A z [b zz ] , must be performed to determine their values. A more compact, “grid” 
type stencil is defined by multiplying Eq. (A. 2) by A zz A yy A xx : 

(A. 3) [ A z z AyyB vr + A zz A xx Byy + Ayy A XX B ZZ \P — A zz AyyA rx F 

where the commutation properties, A yy A. rr = A xx A yy , A ZZ A XX = A XX A ZZ , and A zz A yy = A yy A zz 
have been used. Using tridiagonal compact finite differencing for the matrices in Eq. (A.3) results in an 
explicit 27-point stencil for the LHS and RHS of Eq. (A.3). Because of the symmetry of the second deriva- 
tive, there are only 8 unique coefficients in the RHS and LHS of 27-point stencil which are given in Table 
A.l. Tilt' stencil represented by Eq. (A.3) does not require a separate 1 near boundary scheme for Dirichlet. 
Ixmndary conditions since the stencil is three points wide. 
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Table A.l Coefficients of the LHS 2 /-point explicit stencil for 3D Poisson equation using fourth-order 

compact scheme. 


Location 

Coeff.s of LHS 

Coeff.s of RHS 


-2a(l/h 2 x + l/h 2 y + l/h 2 z ) 

1 


a(l/h 2 x - 2a/ h 2 - 2a/ h{) 

a 

i,J± 1> k 

a(-2a/h 2 + J/fi 2 - 2a/ h 2 ) 

a 

U j,k.± J 

a(-2a/hl ~ 2a/h l + I/h b 

a 

i ± 1, j ± 1, k 

aa(l/h 2 x + l/h 2 y -2a/h 2 z ) 

2 

a 

i ± 7, j y k ± 1 

aa(l/h 2 x ~2a/h 2 + 1/h 2 ) 

2 

a 

i, j±l,k±l 

aa(-2a/h 2 r + l/li 2 v + l/h 2 ) 

2 

a 

i±l, j±l,k±l 

aa 2 (l/h 2 x + 1/h 2 + l/l 2 ,) 

3 

a 


APPENDIX B - Solution procedure for curvilinear grids. The procedure for solution of the 
governing equations on general grids using compact finite differencing is outlined in the this appendix. The 
approach taken is to transform the spatial gradients from physical space to a computational space with 
uniform grid spacing where higher-order, compact finite differencing is used. The velocity components are 
not transformed and are defined in the Cartesian coordinate system. 

First derivatives in the Cartesian coordinate system are expressed in terms of derivatives in the compu- 
tational coordinate system using the chain rule: 


(B.l) 


3<f> 




3<t> 


dx j dx j de m 


Xlic Laplacian operator ill the Cartesian coordinate system is expressed in terms of gradients in the 
computational space: 


(B.2) 


d(p 


dc 


dx-dx- dxjde n 


a 2 <t> 

l,, dE,de I 


+ Til 


r<Kd j> 

3 X j dE n 

a 2 $ 

22 de 2 de 2 


+ m 


a 2 <t> 




3<|> 


> 2 dE,de 2 +m ’dr ] +m2 dr 2 


where: m }] = 




m 


22 


= fixl + yl) 


m 


12 


= - 2j2 [ x t\ + y e y^ » 


./ = 


( x e^ 


yt x r\) 


m l = J ^-y^ m n x tt + ln 22 x r\r\ + "*«*«!,) + x r,(™/^ee + m 2 2 V^ + „)] 

m 2 = J[y t (m n x tt + m 22 x m + m /2 x en ) - x z {m n y zz + m 22 y ^ + m /2 j/ eT) )] 
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Spatial gradients of the computational coordinates appearing in Equations (B.2) are the metrics from a 
direct transformation and are defined in terms of the indirect metrics: 


x t 'n x \ 

y E v$ 

_ 2 e F \ 

Compact finite differencing is used to evaluate the gradients of the Cartesian coordinates in computa- 
tional space that appear on the RHS of Eq. (B.4). 

Equations (B.2) and (B.3) are used to express spatial gradients appearing in the continuity and 
momentum equations. The continuity equation (2) in the computational coordinate system is given by: 


(B.3) 


e .r e v C. 




3e, 3 u : 


(B.4) 


■' = 0 
3.^3^, 


while the momentum equation 

(1) is given by 




8 ", de „, du i 

m d 

f3£„3H, 

(B.5) 

dt + U ->dxj 3e„, 

~ ~dx i de„, h Re jjdxjdE,,, 

l d -‘j de n 


Discretizing Eq (B.G) temporally with the RK scheme and spatially with compact finite differencing 
gives: 


(B.G) 


M + 1 M , M 

u- = u- + t) 




via 


re Hf = - «f(8,.£„,)(6 Ei( «f) + ]^( 7 »ii 8 ee + + ™ /2 8 eii + m A + m 2 8 n )w ' W 


Taking the divergence of Eq. (B.G) gives the Poisson equation for pressure. 


(B.7) 


(”*// 8 ee + m 22 8 nu + TO /2 8 en + 711 ih + m 2 8 u^ 

= ( 8 r. e »») 5 E 


M 


M 


h m + 


b‘ XI+I At 


Explicit stencils for Eq. (B.7) are difficult to find due to the presence of the metric terms which prevent 
commutation of the derivative operators. Solution techniques for Eq. (B.7) are based on treating the pres- 
sure, F W , as well as it second derivatives, 8 ee F W and 8 ni) F W , as unknowns to be solved simultaneously: 

(B.8) ' 4 E E ( 8 EE /yU ) = B tt 
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(B.9) 


^„„(5 riT ,/ yU ) = B P^ 1 
Tin v nn ’ uu 


(B.10) 


(m // 8 ee + m 22 8 rin )/ y 


M 


= ( 6 x, e m) 5 e 


\I u 

Hr + 


M 


M 

b At 




- ( m /2 5 eri + ™/ 6 e + m 2 8 Tj)^ 


where A EE , B (y . .1^ . and are the compact matrices for the second derivative scheme defined in 
Section 2.3.2. The first and mixed derivative terms have been placed on the RHS of Eq. (B.10) and are 
lagged at the previous iteration. Iterative techniques outlined in Sect. 3 are then used to solve the system 
of equations (B.8) - (B.10) for the pressure and its second derivatives. 
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litude. IRI 


Fig 5 




Solution to the 1-D convection equation at t = fOO in wavenumber space for various finite difference 
approximations of the first derivative term f (a) amplitude and (b) phase angle. 
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Fig 8 




Numerical solution of the rotating cone problem after one revolution on a 82 x 82 grid (a) third- 
order upwind scheme, (b) fourth- order compact scheme. Numerical solution is shown to the left , 
exact solution to the right. 
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Fic 9 




Numerical solution of the rotating cone problem after one revolution on a 32 x 32 grid from Orszag 
(1971), (a) second-order Arakawa scheme , (b) fourth-order Arakawa scheme, and (c) spectral meth- 
ods. 
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Solution of the Stuart's problem, (a) Exact solution and (b) numerical solution using fourth-order 
compact scheme on a 1,1 x 41 grid at t — 20jt 
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